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Abstract. We present a pulmonary vessel segmentation algorithm, 
which is fast, fully automatic and robust. It uses a coarse segmentation 
^^ of the airway tree and a left and right lung labeled volume to restrict a 

vessel enhancement filter, based on an offset medialness function, to the 
I— I lungs. We show the application of our algorithm on contrast-enhanced 

^>' CT images, where we derive a clinical parameter to detect pulmonary 

r\ hypertension (PH) in patients. Results on a dataset of 24 patients show 

that quantitative indices derived from the segmentation are applicable to 
distinguish patients with and without PH. Further work-in-progress re- 
sults are shown on the VESSEL12 challenge datasetjj which is composed 
of non-contrast-enhanced scans, where we range in the midfield of partic- 
ipating contestants. 
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1 Introduction 
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_|! Since its introduction in the 1970s, computed tomography (CT) has become 

^^ an important tool in medical imaging. It is the gold standard in the diagno- 

CO sis of a large number of different disease entities ^, and further technological 

^"^ progress has strengthened its diagnostic impact leading to an essential role in 

^ clinical practice. To gain full benefit of the increasing resolution of CT im- 

V^ ages, automatic methods are needed to separate important information from 

rN diagnostically irrelevant ones. 

^ Automatic segmentation and analysis of the vessels inside the lung (pul- 

monary vessels) from CT images is widely used for computer aided diagnosis 
of vascular diseases [17 , non-rigid image registration [13 , and detection of pul- 
monary embolism [18 . Our clinical focus is on the detection of pulmonary 
hypertension (PH), which is a chronic disorder of the pulmonary circulation, 
marked by an elevated vascular resistance and elevated mean pulmonary artery 
pressure (mPAP) [7]. Unlike the systematic circulation, the blood pressure in 
the pulmonary vessels is very difficult to measure. In order to determine pul- 
monary pressure, an invasive right heart catheterisation must be performed [7]. 
By the time of diagnosis, PH has usually progressed to late stage and is not 
reversible any more. By finding a non-invasive way of measuring the pulmonary 
blood pressure, the number of patients awaiting treatment could be significantly 
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decreased. Therefore, an important aim of clinical PH research is the early di- 
agnosis of pulmonary hypertension. 

For determining a measure for pulmonary blood pressure, a segmentation 
of the blood vessels inside the lung is needed. We present a robust algorithm 
that combines lung- and airway-segmentation, together with a sophisticated 
vessel enhancement filter to obtain a proper segmentation of the left and right 
pulmonary vessel trees, even in patients showing severe pathologies. A schematic 
of the anatomy of the human lung is shown in Figure [l] 

The algorithm is fully automatic, computationally efficient and able to han- 
dle large datasets. It was tested on phantom data, the publicly available VES- 
SEL12 challenge dataset, and CT data from 24 patients from a clinical PH 
study. 



Trachea 




Figure 1: Anatomy of the human lung showing the airways, starting with the 
trachea, left and right lung, the pulmonary arteries in red and the pulmonary 
veins in blue (© www.somersetmedicalcenter.com). 



1.1 Related Work 

Several 3D vessel segmentation algorithms have been presented in the literature 
up till now. An overview can be found in |12j. Typical approaches include 
threshold-based algorithms [9 or fuzzy methods [8]. These approaches have in 
common that an intensity model is utilized to detect the vessels. Frangi et al. [6] 
presented a technique based on the eigenvalue analysis of the Hessian matrix. 
This approach was later refined in the popular approaches of [17^ and [10^, who 
also take the eigenvectors of the Hessian matrix into account. In [1] they ap- 
ply vessel-enhancement, junction-enhancement and nodule-enhancement filters, 
based on an eigenvector analysis from the gradient vector field. Most of these 
approaches work well on controls, however, for patients showing pathologies, 
robust extraction of vascular structures is still an open issue, especially in the 
case of pulmonary hypertension where vessel pruning occurs [11^ . 



2 Method 

At the core of our method is a multi-scale vessel enhancement (VE) filter based 
on the Hessian matrix. It is similar to [15 in using the eigenvectors of the 
Hessian matrix to detect candidate voxels inside the vessels, and computing 
an offset-medialness boundary measure perpendicular to the estimated vessel 
direction [10]. The VE response (i.e. medialness) is limited to the right and left 



OAGM/AAPR Workshop 2013 IjarXiv: 1304.187^ 



lung, which is derived from an intensity-based lung segmentation and a coarse 
airway tree segmentation with a labelling of left and right main bronchi. After 
centerline detection from the VE response, a coarse radius estimation of the 
vessel is obtained using a spherical ray-casting approach. The final segmentation 
is the output of a globally optimal geodesic active contour model based on total 
variation and a data term derived from the medialness [16] . Figure [2] shows the 
fiowchart of our automatic vessel segmentation approach. 




Figure 2: Flowchart of the vessel segmentation algorithm, (a) input CT-image, 
(b) lung- and airway segmentation, (c) vessel enhancement filter restricted to 
the lung (yellow: high vessel probability), (d) centerline reconnection, (e) vessel 
segmentation 



2.1 Lung and airway segmentation 

A prerequisite for our vascular tree extraction is a segmentation of left and right 
lungs, react ively, to restrict the reconnection of the vessel centerlines. A coarse 
airway segmentation and labelling of the main bronchi initiates this process. 
The labeled airway tree is subsequently used to label a coarse, threshold based, 
lung segmentation. Segmenting the airways also helps in removing false positives 



of the vessel enhancement filter (see Section 2.2), since the intensity contrast 



of the airway border and blood vessels is very low, thus leading to incorrect 
detection of blood vessels at the airway walls. 

We automatically detect the airway on the top-most slice of the contrast- 
enhanced volume, which is a dark circle surrounded by high- intensity tissue, to 
get a seed point for an iterative 3D region growing algorithm. For the region 
growing two thresholds are defined: thmin = ^(^s) — 1 HU and thmax = H^s) + 
1 HU, where / is the CT image, Xg is the seed point and HU denotes a Hounsfield 
Unit. All A^6 connected voxels which fulfill thmin < /(x) < thmax are added to 
the segmentation. Then the thresholds are updated {thmin = thmin — 1 HU and 
thmax = thmax + 1 HU) and region growing is restarted with the new segmen- 
tation as seed. Convergence of this iterative procedure combines two stopping 
criteria. We check whether the number of total voxels, or the number of edge 
voxels of the current segmentation is three times larger than the average of all 
of the previous voxel counts. In this case a leakage is detected (see Figure Isb), 
and we produce the final airway tree segmentation with the restricted thresh- 
old range from the previous iteration. Figures [3^ shows an example of airway 
segmentation obtained from one patient in the clinical PH study. 

After coarse airway segmentation, we perform a left and right lung segmen- 
tation to identify a region of interest for later vessel detection. A coarse lung 
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segmentation is obtained using Otsu's optimal thresholding method [U]. With 
a connected component analysis, the lung is selected and a 3D hole filling is 
applied to include vascular structures. The two lungs always merge through the 
airways, but in some datasets the border between right and left lung is hardly 
visible, resulting in connected lungs. The airway segmentation is used to sepa- 
rate the coarse lung segmentation. Using a graph representation of the skeleton 
from the airway segmentation, we detect the carina (where the trachea splits 
into the left and right main bronchi), and assign different labels to the trachea, 
right and left bronchi. To label the voxels in the coarse lung segmentation, we 
calculate shortest paths to the labeled airway tree, thus, splitting it into left 
and right lung. As a cost function Ic for the shortest path algorithm, we use 
the gradient magnitude of the CT image |V/| and the coarse lung segmentation 
without the airways, where we give a larger weight to the gradients to prevent 
connecting via short cuts from touching left and right lungs. We found a frac- 
tion of I for the gradient weight, and a weight of ^ for the lung segmentation 
without the airways to perform well in the labeling of our datasets. 

As a final step, to remove holes caused by vessels and other high intensity 
structures inside the lung, the airways are removed from the lung segmentation, 
and a morphological closing operation is applied several times at each lung 
separately. We use a six-neighbourhood star-shaped structuring element and 
10 closing operations. These parameters remain constant for all datasets. This 
ensures that the lung segmentations contain the whole lung. The different steps 
of the lung segmentation can be seen in Figure |4] 




(a) 




Figure 3: (a) representative result of the automatic airway segmentation from 
an example patient, (b) leaked airway segmentation 




(a) 



(b) 



(c) 



Figure 4: (a) one 2D slice of an example CT image, (b) coarse lung segmentation 
after separation, (c) refined lung segmentation, separate left (white) and right 
(grey) lung 
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2.2 Vessel enhancement 

We enhance all vessel-like structures using a modified version of the vessel en- 
hancement (VE) filter proposed by Pock et al. [15 . We extend the VE filter by 
using multiple radii in addition to a multi-scale pyramid with a down-sampling 
factor of 1.7. With these parameters, the experiments showed improved dis- 
crimination between different vessel radii compared to the widespread use of 
Gaussian pyramids with a down-sampling factor of two. Further, to be more 
robust against noise and to get rotation invariant derivatives, we compute our 
gradients for the boundary measure and for the Hessian matrix using the filter 



kernels of [5 . The airway- and lung-segmentations from Section 2.1 are used to 
restrict the vessel enhancement output to the left and right lungs without the 
airways, respectively. 

To get the vessel enhancement filter response, we calculate the eigenvalues ei, 
62, and 63 (ordered such that |ei| > |e2| > leal) and the associated eigenvectors 
vi, V2 and V3 of the Hessian matrix H^(x) at each scale a. To sort out for 
bright tubular structures on dark background, we check that ei < and 62 < 
holds. In points that fulfill this condition, the smallest eigenvector V3 gives 
an estimation for the vessel direction. Perpendicular to the vessel direction, 
in the cross section plane of the tube given by the eigenvectors vi and V2, 
we evaluate boundary information along circles of different radii r. We define 
the boundary gradient B(x, a) = (tV/^(x), with /^(x) being the CT image 
convolved with a Gaussian kernel with variance a. An initial response is given 
by the median of the boundary contributions bi^a- = |B(x + tVq,., ctJVq,. |, with 
i = 1 . . . [27rr + Ij , which we denote as Rq . A problem of i^o^(x, r, a) is, that it 
also produces responses at isolated edges. To avoid this, a measure of symmetry 
is introduced: 

6(x, r) = 1 = — 



where 5(x, r) is the median absolute deviation of the boundary samples and 
b is the median. The final boundary response is computed as: 

i?+(x,r,cr) = i?J(x,r)6'(x,r)i 

To suppress responses at the border of vessels, the gradient magnitude at the 
center of the vessel is combined with the offset medialness R^ from above: 

i?(x, r, a) = max {i?+(x, r) - cr| V/^(x)|, O} 

The final vesselness response 

Rmuitii'x.) = max {i?(x, r, a)} 

a,r 

is the maximum response from all different scales a and radii r. We found 4 
scales and radii r varying from 1 to 2 pixels, with an increase of 0.3 pixels to 
have the best performance. 

2.3 Centerline extraction 

In a non-maximum suppression step inspired by [?, at each position x with 
a medialness i^(x) > thmim we sample 8 points on a plane perpendicular to 
the estimated vessel direction. If the medialness on any of those 8 points is 
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larger than at the current position x, the VE response at x is set to zero. 
This leads to disconnected vessel center line fragments due to branching points, 
where the tubularity assumption fails. Next, small centerline fragments (less 
than 5 A/'26-connected voxels) are removed, and all maxima lying on the airway 
border are cleared. To reconnect the centerline fragments, we apply a Dijkstra- 
like shortest path algorithm. In each lung separately, we connect all centerline 
candidate points to the center of the heart. As a cost function, we combine 
the medialness with the gradient magnitudes of the CT image. The separate 
processing of right and left lung ensures that wrong connections through the 
mediastinum are avoided. The connected trees of the right and left lung form 
the final vessel tree. 



2.4 Vessel segmentation 

At each centerline voxel, 48 points lying on a sphere are sampled from the 
CT image and the grey values are summed up. This is done for spheres with 
different diameters. As soon as the normalized sum of grey values drops under 
0.6 (empirically found), the radius of the vessel has been found. This gives a 
more accurate radius estimation compared to the filter radii directly from the 
VE filter response. This coarse segmentation is then used as an input for a total 
variation based segmentation [IGJ, which gives the final segmentation. 



3 Experimental results 

For testing our algorithm we used the publicly available VESSEL 12 challenge 
dataset as well as 24 datasets from patients who underwent contrast enhanced 
CT as a part of a clinical PH study at the Ludwig Boltzmann Institute for Lung 
Vascular Research, Graz. The median size of the CT volumes is512x512x 426 
pixel. Our hardware consists of an Intel CORE I7-2600K 3.40GHZ with 16 GB 
RAM and a Geforce 580 GTX with 3 GB RAM. 

All of our filters are implemented in CUD A, which is a parallel program- 
ming model from NVIDIA, that generates hardware accelerated instructions for 
NVIDIA graphics processing units (CPUs). Using CUDA we can significantly 
improve the runtime of our 3D image processing algorithms, due to the high de- 
gree of parallelization. Therefore, we can directly work on the full resolution CT 
data. The runtime of the CUDA implementation of the whole algorithm pipeline 
ranges from 5 to 10 minutes, whereas a CPU-only implementation needs more 
than an hour. Limited memory on our GPU restricts the possible size of datasets 
in the different stages of our pipeline, so we decompose large (i.e. nearly 512^) 
CT images into overlapping sub- volumes, which are processed sequentially, with 
each sub- volume benefiting from the CUDA based parallelization. 

3.1 Airway and lung segmentation 

On 24 patients we were able to extract the airway tree to at least four gen- 
erations, which is sufficient for our further processing. Due to GPU based 
parallelization, one airway segmentation takes on average 69 seconds. Figure [3^ 
show a representative result. 

The lung segmentations from the 24 patients included all lung tissue and 
pulmonary vessels, which has been verified through visual inspection (Figure [5|. 
Processing time for the lung segmentation was on average 132 seconds. 
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(c) 



Figure 5: (a) 3D rendering of an example lung segmentation, (b) coronal and 
(c) axial view; green=left lung, blue=riglit lung, red=tracliea 

3.2 Phantom data set 

We used a phantom of a liver vessel tree, depicted in Figure |6|i, to check the 
performance of the algorithm and validate its robustness against noise. We 
successively added Gaussian noise with increasing variance to the phantom data, 
and calculated the Jaccard index of the ground-truth segmentation with the 
obtained segmentation. The curve in Figure [6J3 shows how the Jaccard index 
changes if Gaussian noise with increasing variance is added to the phantom 
dataset. As long as the variance of the noise is below 40 Hounsfield Units (HU), 
the performance lies above 93% segmentation overlap. 





100 120 



(a) 



(b) 



Figure 6: (a) 3D rendering of liver vessel phantom, (b) Jaccard index over 
variance of Gaussian noise 



3.3 VESSEL12 challenge 

We applied our algorithm to 20 datasets made available through the VESSEL 12 
challengq^ Table [I] shows the results. The algorithm performs very well in 
terms of specificity, however improvement in the sensitivity is still necessary. 
This is because the algorithm is optimized for contrast enhanced CT images 
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and for finding vessels even in noisy datasets. Thus vessels smaller than 2mm 
in diameter are misclassified as noise and are not included in the segmentation. 
All results from the other participating groups can be found on the official 
VESSEL12 challenge website. 



Table 1: VESSEL 12 challenge results; our team in comparison with the currently 
best performing team (LKEB China) 







LKEB China 




Our Team 




Az 


Specificity 


Sensitivity 


Az 


Specificity 


Sensitivity 


Small vessels/Non-vessels 


0.979 


0.935 


0.956 


0.781 


0.971 


0.585 


Medium vessels/Non-vessels 


0.982 


0.960 


0.957 


0.848 


0.971 


0.716 


Large vessels/Non-vessels 


0.991 


0.983 


0.977 


0.951 


0.971 


0.909 


Vessels/ Airway-walls 


0.940 


0.832 


0.953 


0.851 


0.928 


0.742 


Vessels/Dense abnormalities 


0.683 


0.460 


0.953 


0.712 


0.572 


0.742 


Vessels/Mucus-fiUed bronchi 


0.513 


0.204 


0.953 


0.677 


0.439 


0.742 



Az: Area under the Receiver operating characteristic curve; Specificity: number of true 

negatives / (number of true negatives + number of false positives); Sensitivity: number of 

true positives / (number of true positives + number of false negatives) 



3.4 Contrast enhanced CT images from our clinical PH study 

Our clinical application is the detection of pulmonary hypertension (PH), which 
is a chronic disorder of the pulmonary circulation, marked by an elevated vas- 
cular resistance and elevated mean pulmonary artery pressure (mPAP). Our 
hypothesis is, that the pulmonary vascular tree shows quantifiable differences 
between patients with and without PH. One quantifiable property of the ves- 
sels is their tortuosity, which is a readout of twistedness [2| . The most common 
metric of vascular tortuosity is the distance metric (DM), which provides a ratio 
of the actual vessel length to the Euclidean distance between its endpoints [3j . 
To determine the tortuosity, the lung vessel centerlines and branching points 
are extracted. The DM is calculated and compared with the patient's clinical 
parameters. 

We found a correlation between DM and mPAP of p = 0.60 (Spearman cor- 
relation coefficient, p<0.01). There was a significant difference between the DM 
of patients with and without PH (Table [2J p<0.05), thus enabling to discrim- 
inate the two groups on our dataset of 24 patients. Two representative vessel 
segmentation results of the PH study datasets are shown in Figure [7| 



Table 2: Distance metric 



Readout 



All patients (n=24) 



No PH (n=6) 



PH (n=18) 



Distance metric 



1.207 ±0.017 
(1.184-1.248) 



1.193±0.011 
(1.184-1.212) 



1.211 ±0.017* 
(1.189-1.248) 



Data is presented as meanzbstandard deviation (range). The significance was tested with 
Students t-test; * p < 0.05 as compared to No PH group 
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(b) 

Figure 7: Representative results of vessel segmentations showing a patient with 
(a) and one without PH (b). No visual differences between the vessel structure 
can be seen. 

4 Conclusion 

We have presented a segmentation approach for vascular structures from contrast- 
enhanced CT images using a multi-scale vessel enhancement filter and using 
information from a lung- and airway-segmentation. We achieved very good seg- 
mentation results on our 24 patients from a clinical PH study. We also tested the 
algorithm on non-contrast-enhanced data from the VESSEL 12 challenge, where 
we occupy a midfield position among all participating teams. We see room 
for improvement in the case of small vessels. Reasons for this performance of 
the algorithm are the optimization for the contrast-enhanced setup and the use 
of isotropic CT scans, which is not the case in the VESSEL12 datasets. Due 
to a parallelized CUDA implementation, our whole vessel tree segmentation 
and centerline extraction shows a run-time of at most 10 minutes for large CT 
datasets, without the need for computing on reduced resolutions, thus enabling 
the potential use in daily clinical routine. 

As an important outcome of our work, we showed that tortuosity is correlated 
with mean pulmonary artery pressure, and our vessel segmentation algorithm 
can detect the presence of PH. One of the limitations of this study is the small 
number of patients, which allows only preliminary conclusions. A large scale 
prospective study to determine the true benefits and constraints of this method 
is currently in planning. Further, due to the radiation exposure one cannot 
test the repeatability of the method. This would be necessary to determine its 
ability for use in disease monitoring and follow-up examinations. 
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